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The problem of the estimation of distribution parameters in the case of experi- 
mentally limited information is discussed. As an example the determination of the 
total number of muons in Extensive Air Showers registered in the KASCADE exper- 
iment is studied in details. Some methods based on other than standard maximum- 
likelihood approach are examined. The advantages of the new methods are shown. 
The comparison with the Artificial Neural Network approach is also given. 
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^ r- 1 1 Introduction 

a: 
a: 

The experimental study of a particular distribution shape is often related 
to the extraction of an information from the statistically scant data. Due to 
the physics of a phenomenon under study or to the experimental situation the 
distribution we want to distinguish is only sampled by the very limited number 
of measured occurrences. The very good example of such situation takes place 
in cosmic ray ground level physics. The high energy particle cascade initiated 
by the particle of energy of PeV or more reach the sea level as a huge number 
of elementary particles. More that 95% of them are electromagnetic particles 
(e~, e + or high energy gamma quanta) there are also very few hadrons and few 
percents muons. All these particles are distributed in the plane perpendicular 
to the shower axis. Shapes of these distributions in each Extensive Air Showers 
(EAS) contain an information about two main points of interest in cosmic ray 
physics: the nature of the primary particle and the character of the very high 
energy interaction. Thus, to extract that information the distributions have to 
be regain first. The electrons (positrons and gammas) are relatively abundant 
so the determination of their distribution makes rather modest problem. The 
muons, as it will be shown below, gives an opportunity to look for a non- 
standard solutions of the distribution parameter estimation problem with a 
poor statistics. 
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2 Basic definitions 



Let the f(x; {p m }) will be the "inclusive" distribution of some random vari- 
able x. The form of / is preadjusted from some theoretical assumptions or 
experimental experience. There is a set of parameters {p m } (including may 
be also a normalization constant - the / need not to be normalized to unity) 
to be estimated. The methods commonly used for that purposes are based on 
minimization with respect to {p m } some distance measure between the n par- 
ticular random realizations of variable x and the predictions for given values 
of {pm}- (The value of n can be also a random variable, as it will be seen in 
the case of muons in EAS.) The standard procedure is to build the histogram 
from the {xf, i = 1, . . . , n} and use one of handbook methods to fit the his- 
togram using values obtained from f(x; {p m }) by respective integration. The 
most popular method is to use the y 2 measure. 

The problem can be solved satisfactory in all the cases where the method is 
applicable. In general it can be said than this happens if the number n of 
used random realizations of function / is sufficiently high. The meaning of 
"sufficiently" is not well define. All depends on the way of binning. For too 
many (with respect to the value of n) too small bins the contents of some bins 
can be not enough to satisfy the x 2 applicability conditions. It is quite clear 
that if we have very few entries in a particular bin the fluctuations of the bin 
contents could not be treated as gaussian and the simple probabilistic meaning 
of x 2 is no longer valid. However the formula for calculation the distance: 

_ (Nei-Nci) 2 

U ^ £ jy" > C 1 ) 

histogram * 

where Nei is the content of z-th histogram bin and Nci is its expected value 
for a given set of parameters, can be still treated as a measure of a distance 
between the histogramed data and the prediction for given {p m } remember- 
ing that the probability distribution of the variable u is not the x 2 one - m 
that sense it can be used in minimization procedures to fit some unknown 
parameters for any bin contents. 

The inconvenience of statistical interpretation of the distance defined in Eq. 
([]]) can be avoid if we use the exact maximum likelihood method. The measure 
used is: 

u := £ [ — In ( p (-/Ve^; iVc*) ) ] , (2) 

histogram 

where p(Nef, Nci) is the probability that i-th bin of the histogram contents 
is iVej while the value expected for the given set of parameters {p m } is iVcj. 

In that paper the other possibilities of the distance measure between his- 
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togramed data and given values of {p m } prediction are also investigated. The 
first proposition is the measure similar to the euclidean metric in multidimen- 
sional space: 



u 



the other is: 



/ ]T (Ne t - Nc t ) 2 , (3) 

histogram 



u := ]T \Ne t -N Cl \. (4) 

histogram 

During the binning procedure some amount of information is lost. The problem 
arises in particular when the normalization constant is one of the parameter to 
be estimated and the expected number of n is small. The methods which use 
the information about each single x (or from histograms with very small bin 
size, what is analogous to some extend) establish qualitatively quite different 
and important (as it will be seen) class of measures. Among them the best 
known are those derived using the concept of random walk (which is, in fact, 
a basis of the well known Kolmogorov-Smirnoff test). From the set of {xi} 
something like the cumulative distribution can be made: 

k 

De j = IZ 9 ^-^), (5) 

i=i 

where 

„ , , J for x < , s 

0(X) = (l for x>0 (6) 

and j can enumerate both: single values of x or histogram channel contents. 

In the same way the expected "cumulative distribution" (normalized with 
respect to the n) of f(x; {p m }) obtained for given values of parameters can be 
defined: 

DC 1 <$>(Xj ~ x)f(xAPrn})dx, (7) 

j 

where the integration is taken over the whole x domain covered experimentally. 
Using both "cumulative distributions" the measure can be defined as: 

u := max ( \Dej — Dcj\) (8) 

Here again some modifications can be made. Two examples will be discussed 
in that paper: 

u := max (Dej — Dcj, 0) + max (Dcj — Dej, 0) (9) 
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and 



u 



max (Dej — Dcj , 0)^ + ^max (Dcj — Dej, 0) 



(10) 



3 Applications to the muons in KASCADE experiment 



For the statistical methods examination it is enough to test the case of ideal 
experiment. Particle detectors are distributed over some area and each detector 
responses giving the number of muons passed through its surface. In the real 
situation there is some spread of the detector signal so the results obtained in 
the present work can be considered as a most optimistic approximation of the 
reality. 

The geometry of our ideal experiment studied in that paper is exactly the ge- 
ometry of the KASCADE experiment (Ref. Q). This experiment is at present 
about starting to collecting data and the expected quality of these data gives 
a hope to make a great improvement in our knowledge of the nature of EAS. 
Among the others there are 192 muon detectors of area 3.2 m 2 distributed over 
the area of 200 m x 200 m. The KASCADE experiment is dedicated to EAS 
induced by primary CR particle of energies from about of 10 15 eV. In that 
paper the low energy threshold of the experiment due to lack of information 
about muon distribution from the array detectors is investigated. 

To start any ideal experiment a source of ideal inputs is needed. The perfor- 
mance of the data evaluation can be tested only when there is a set of presumed 
inputs as they are expected in real case and the respective set of true answers 
which the experiment should give. In our case the Monte-Carlo shower simu- 
lation code CORSIKA version 4.112 was used (Ref. [Q]). That program realize 
the complete simulation of particle passage through the atmosphere producing 
the large output file of all particles on the level of observation. For that paper 
the use of particular code is not very important. The statistical behaviour of 
an ideal experiment does not, to some extend, depend of the physical exactness 
of the Monte-Carlo calculation of EAS development. For the same reasons in 
the present analysis only vertical showers are used. 

The Monte-Carlo program in principle do not give the muon lateral distribu- 
tion. Some finite number of muons is nevertheless distributed somehow over 
the experiment surface. Due to the random character of physical processes in 
the shower development the concept of muon lateral distribution in each indi- 
vidual EAS can be introduced. It can be define as a probability distribution 
sampled by each muon in Monte-Carlo simulated shower. All the information 
about that distribution can be obtained, by definition, only from the finite 
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number of all muons in the particular shower. The accuracy of such definition 
for our purposes is quite enough (the surface covered by muon detectors allows 
us to measure only the fraction of percent of the muons in EAS). 

The detail shape of the muon lateral distribution is of course unknown, but 
some analytical formulae could be found in the literature. The particular choice 
is again not very important for the present study. We used the one which 
reasonably well describes the Monte-Carlo outputs: 



where r is the distance from the shower core, which position was assumed to be 
known for our purposes (from electromagnetic EAS component measurement). 
The normalization of the above distribution was done in the way that the p 
can be treated as a muon density at a given distance and is the total muon 
number in the shower. There are four parameter in Eq. ([TT| ) to be adjusted 
(Np, Ro, n and m). It has been checked that the value of m can be fixed without 
loosing the quality of the ideal inputs description. The value of m = 1.7 is used 
hereafter. This was possible mainly because of the fact that the distances from 
the shower core to detectors were limited to about < 300 m (the simulated 
showers were uniformly distributed within about 50 m form the center of the 
KASCADE array). 

For each simulated shower the formula in Eq. ([TT]) was fitted with the help of 
standard methods to all muons in the distance range from 10 m to 300 m. The 
exactness of such fits is given below. The obtained parameters N^, R and n 
will be hereafter called true parameters of the particular simulated shower. 

In the next step the muons from the simulated shower are sampled by the net of 
our ideal muon detectors producing the output signal of the ideal experiment. 
Detectors can be treated as a relatively small bins in the histogram with 
respect to the distance from the shower core. 

The general problem is how the information about true values of the individual 
shower muon distribution parameters (JV M , Rq, n) can be derived from the 192 
histogram channel contents. 

First method tested was the standard \ 2 method. The histogram (detector 
responses) was rebind. Finally, after some tests, the constant bin limits have 
been chosen: (10, 20, 30, 50, 70, 100, 130, 170 and 200 m). The sum in Eq. 
([]]) was performed only for the bins with the contents not less than 5. That 
method will be hereafter called Hxl. 

Without assumption about the minimal bin contents the Eq. (JJ) was also used. 
This method will be called Hx%- 




(11) 
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Next tested possibility is associated with the Eq. (|2|) The assumption that all 
the muons are uncorrelated (for one particular EAS) leads to the poissonian 
distribution of the fluctuations in each histogram bin: 

Yd 1 — 

p(n; n) = — e~ n . (12) 
n! 

Using that directly with Eq. (fj) the method called HL is defined. 

The equations Eqs. (||) and (§) applied to the histogramed detector outputs 
leads to the methods labeled HE, HA respectively. 

The Kolmogorov-Smirnoff like method was also checked for the rebind his- 
togram. The Eq. (H) where j indexes 10 big bins lead to the method called 
HK. 

The main aim of that work is to look for the best estimation method. It 
was foreseen that using directly each detector information the results should 
be improved. First the standard measure (Eq. (|l])) was checked. ( With Ne^ 
replaced in that case by the i-th detector response and Nci by its expected 
value. ) 

The assumption of minimum hits in particular detector was of course rejected. 
That method will be called D%. 

With analogy to the method introduced to the big bin histogram the measures 
defined by Eqs. (|j) and (^) were used leading to the methods denoted by DE 
and DA respectively: 



u 



{rii-s p(ri)) 2 , (13) 

detectors 



u := \ n i ~ s P( r i)\ > ( 14 ) 

detectors 

where s is a detector area, its distance from the shower core and rij detector 
response. 

The maximum likelihood methods are the standard tools in the estimation 
theory. In the present work that possibility was investigated too. In general, 
it is based on the knowledge of the shape of fluctuations with respect to the 
mean of the muon number passing the given detector surface at fixed distance 
form the shower core. The detail study of the Monte-Carlo simulated showers 
allows us to state that they are wider than poissonian (Ref. ||). They can be 
well approximated by the negative binomial distribution: 



p(n; n,j) 






1/7 


n 
7 


1/7 


') 


[ 1 + 1/7. 




[l + l/7 J 



(15) 
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with 7 = 1.0. With that assumption the maximum likelihood measure can be 
defined (and than minimize): 

u := ]T [-ln(p(rn; s p(n), 1.0))]. (16) 

detectors 

That method will be hereafter called DPI. 

To see if the particular negative binomial shape of the fluctuations produce an 
important effect the measure defined by Eq. (|2|) was used with the poissonian 
fluctuations (Eq. (|P2|)). The results of such a calculations is denoted by DP2. 

The Eq. (§) is of particular interest for the individual detector data. The 
minimization of that measure leads to the method called hereafter DD. The 
modifications described in Eqs. (H) and (|ToD were also tested. First named 
DDI and the second DD2. 

The last method results of which are presented in this paper is based on com- 
pletely different way of solving the estimation problem. This is the Artificial 
Neural Network (ANN) method. The much more detailed discussion of it will 
be given elsewhere (Ref. ||). 

For the present analysis the neural network contains 192 input nodes which one 
steering with an integer response of one detector in the array. There are two 
hidden layers with successively decreasing numbers of neurons and finally one 
output node. The output signal of the network is related to one of the param- 
eters (Nfj) of the muon lateral distribution (Eq. [TT| ). The number of weights 
of the network to be adjusted is very large so for the network training about 
hundreds of thousands showers were used. A special technic was developed to 
achieve this. Details are not very important for the present work. Finally we 
obtained the network trained with the proton induced vertical showers in the 
energy range of interest. The ANN method results presented here are given 
just to compare them with the others and should be treated as, more or less, 
preliminary. 



4 Results concerning the total number of muon estimation in the 
KASCADE experiment 

The parameter of the great importance for the EAS study is a total muon 
number in the shower N^. 

To perform the minimization of the respective u measure the standard MI- 
NUIT from the CERN library was used (Ref. ||). It is obvious that in some 
cases due to fluctuations in the input data the minimum of some distance mea- 



7 



a. 



— ANN training sample 
-- 10 14 eV 1 



10 15 eV 
10 16 eV 



1 .. 



10 



1 



N, 



1 



Fig. 1. Distribution of true Nn values for showers of different primary proton en- 
ergies 10 14 , 10 15 and 10 16 eV. The thin dashed line shows the muon shower sizes 
distribution of the artificial neural network training sample. 

sure can not be reached within reasonable limits of the parameters in Eq. flnp. 
The limits used were 10 -j- 800 m for Rq, 0.05 -v- 1.99 for n. The minimization 
was started with two parameters fixed: Rq at 300 m and n at 1.43 (the values 
about the mean for 10 15 eV) and the minimization was performed with only 
one parameter to minimize (iV M ). Then the n parameter was liberated and the 
minimization stared again from the just found values. If possible, next, the Rq 
parameter was fitted also. 

Simulated showers which produced all in our ideal array response matrix 
are of course lost. 

All methods were examined for three different primary cosmic ray particle 
energies. As it has been said the KASCADE experiment is designed to study 
the cosmic ray primary spectrum in the very interesting and unexplored up 
to now region starting at about 10 15 eV. Just below that energy its the end of 
the area covered by the recent top of the atmosphere experiments data. The 
correspondence of the direct (balloon-borne) and indirect (EAS) methods is 
one of the important point of that experiment so the determination of EAS 
parameters at the lower part of the shower size spectrum is very important. 

The distributions of true total number of muons in showers of the energies of 
interest are given in Fig. [I]. 

To better recognize the problem the typical muon lateral density distributions 
of the showers initiated by the primary vertical proton of the energies from 
10 14 eV to 10 16 eV are given in Fig. |^. The true muon distributions (fits of 
form given by Eq. (|TT|) to the histograms) are also presented. 

The spread of the points in the Fig. |] is a result of the physical fluctuations 
in the shower development. To make this figure all the muons in the showers 
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Fig. 2. Examples of the muon lateral distributions for individual proton showers of 
different energies (histograms). The curves represents the true muon lateral distri- 
butions of that particular showers. 



10 14 eV 




10 15 eV 




10 16 eV 




Fig. 3. Examples of the muon response array for showers in Fig. 2. 

were used. The real problem arise when distributions like that are sampled 
with the net of detectors which covered only few % of the shower front area. 
The particular detector array responses for the showers which muon lateral 
distributions were given in Fig. |^ are presented in Fig. [| 



It can be seen that the expected numbers of muons seen in detectors is rather 
small as well as the fraction of fired detectors both for 10 14 eV and even 
for 10 15 eV energies. The situation for the energy 10 16 eV looks better. It is 
interesting that even for the worst case some of discussed methods could give 
the reasonable estimation of the total muon contents. 



The value of obtained for each shower in the ideal experiment using dif- 
ferent methods described above has to be compared with its true value. The 
collection of figures presented below shows the spread of the calculated values 
about the true ones. The number of showers used for each test was about 
3000 and the histograms were normalized in such a way that the area below 
each one is related to the fraction of cases were given method was used. The 
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Fig. 4. The spread of the estimated value for each method discussed in the paper 
for the vertical 10 16 eV proton induced showers. 



shapes of the histograms give a self-explainable answer to the question about 
the exactness of each of the methods. The labels gives the references to the 
particular method for which the histograms were obtained. 

In Fig. |] there are shown results for 10 16 eV. Figs. |5| and ^| present quality 
of the methods for showers with energies smaller of one and two orders of 
magnitude respectively. 



The more quantitative results for different methods comparison are given in 
the Tables. In the 2 nd column mean deviations of logarithm of the calculated 
Nfj, from its true value are shown. In the 3 rd column the accuracy of each 
method defined as a a of the (log 10 Nc^/Ne^) distribution (as they are pre- 
sented in Figs. |]-|6|). The applicability of each method defined as a fraction 
of showers for which particular method was successfully used is shown in the 
4 th columns. In the 5 th column in the Table ^| the fraction of events for which 
at least two parameter fit was possible is given. 
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Fig. 5. The spread of the estimated value for the 10 15 eV proton induced showers. 

For the neural network the total muon number iV M was only one parameter 
to study and a value of it was created for each shower except "all 0" showers 
which were not used in the present analysis. 



5 Discussion of the results concerning the total number of muon 
estimation in the KASCADE experiment 

Staring from the highest energy used for present examinations, 10 16 eV, some 
interesting statements can be given about the efficiency of different methods 
for relatively muon rich showers (N^ ~ 10 5 ) where in every case at least about 
50 detectors were hit. As to the methods based on a histogramed data it is 
seen that all of them lead to very similar results. The reason for that is clear. 
The statistical significance of the information stored in the big bin histogram 
is so large that the particular differences between different measures used in 
minimization procedures has no great effect on the results. The situation is a 
little different when using the raw detector data. The biases are seen for two 
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Fig. 6. The spread of the estimated N„ value for the 10 14 eV proton induced showers. 
The thin lines show results for events for which at least two parameter fit has been 
successfully applied. 

measures: D% and DA. Among the others it is clear that the two using the 
" cumulative distribution" metrics leads to the smallest spread of the calculated 
value with respect to the true one. 

Comparison of the best minimization methods (e.g. DD or DDI) with the 
artificial neural network approach leads to the conclusion that no significant 
difference is seen. This can be treated as an evidence for the assumption 
that the exactness of the value estimation reaches its limit allowed by the 
physical fluctuations of the shower development. If two such different ways 
of the data evaluation shows as good agreement that with the high level of 
confidence one can postulate that it is a real limit, and there is no other way to 
get more information (about the parameter under the study) from that data. 

Going down with the energy to 10 15 eV, it is clearly seen in Fig. |5] and the 
Tables that the quality of all methods decreases. When the number of regis- 
tered muons is really a few per event some methods fail in general but, what 
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Table 1 

Detailed results for 10 16 eV proton shower sample, 
is surprising, some are still useful. 

The differences of using the histogramed data measures appear. Some of them 
become unapplicable in a fraction of events for which the others are still useful. 
The larger differences arise for the raw data. The method DA for such showers 
can not be used any more for showers of that energy. The bias seen for 
method become so great that is usefulness has to be also questionable. It is 
very interesting to note the difference between the two "maximum likelihood" 
method DPI and DP2. The assumption about the poissonian fluctuations of 
the number of muons in single detector leads to enlargement of the error of 
the DP2 method in comparison with DPI. The reason for that is in the events 
in which the relatively big fluctuation on one detector appears. Its probability 
calculated from the poissonian distribution is much smaller that it is in real 
shower situation approximated by the negative binomial density fluctuation 
distribution. The significance in the minimization procedure of that particular 
detector increases when the mean muon density at a detector site is small. 

The shapes of the most of thick histograms in Fig. |^ is not gaussian. This 
is simple a result of the details of the minimization technic used. The initial 
values of the Ro and n parameters were taken to be the average values for 10 15 
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Table 2 

Detailed results for 10 15 eV proton shower sample. 

eV showers. Thus, for a very small showers when the information about the 
muon distribution is very limited those parameters sometimes could not be 
adjusted further. Because the small showers are in average distributed wider 
in shower front plane the total muon number obtained from the fit of only one 
parameter has to be shifted to the smaller values. In the Fig. || there are also 
shown by the thin line histograms the results for the showers for which at lest 
one of the distribution shape parameters could be adjusted. Those histograms 
are much more gaussian. 

Among the minimization methods which works well for those small showers 
it is hard to chose the one significantly better than the others specially when 
trying to compare events with at lest two parameters fitted. 

The method based on the "cumulative distribution" measures (both: using the 
histogramed data - HK, and that using raw data - DD, DDI ) could be used 
for at least two parameter fit in the largest fraction of events. The width of 
the (log 10 Nc^/Ne^) distribution is the smallest for the DPI method. ( Tab. 

I) 

Some advantages of ANN appear much clearly. However the ANN method 
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U.OO 
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U.o41 
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DDI 


-0.081 


0.311 


0.89 


0.309 


DD2 


-0.175 


0.344 


0.89 


0.260 


ANN 


0.331 


0.112 


0.89 





Table 3 

Detailed results for 10 14 eV proton shower sample. 

for that small showers introduces systematical bias toward higher values. 
The reason of that is clear. In Fig. [1] there is presented true shower muon size 
spectrum for training sample. The network was trained only with the showers 
for which at least three detectors were hit while most of the showers of 10 14 eV 
fire not more than five detectors. The value of the bias is correlated with the 
number of fired detectors (Ref. |4j]) so it can be removed "by hand" from the 
ANN output. There is another way of removing that bias. It is expected that 
if there will be used some real physical trigger for the ANN shower training 
sample and the same trigger will be used later in the ideal experiment tests 
the bias will be automatically removed. In such a case even the all zero events 
should produce a reasonable (to some extend) ANN answer. 



6 Conclusions 

Final conclusions of the present work can be formed in a few statements: 

- For large showers the most promising are the methods based on the min- 
imization procedure of the "cumulative distributions" inspired measures 
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using individual detector data (DD and DDI) . 

- For small showers some of the methods using the "cumulative distributions" 
are slightly better than the others. 

- The knowledge of the fluctuation on each particular detectors allows one 
to use right the maximum likelihood method. That knowledge improves to 
some extend the accuracy of the parameter estimation. 

- For the studied in the present paper ideal experiment of the geometry of 
KASCADE (its array muon part) the best developed methods allows one to 
go down with the muon sizes of analyzed showers to about few thousands 
per shower (below 10 15 eV). 

- The artificial neural network approach gives a hope that the limit can be 
shifted down to 10 14 eV. However that methods is strongly effected by the 
exactness of the simulations used for the network training. 
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